Accurate emulation of steady-state and dynamic performances of PEM fuel cells using simplified models

The current effort addresses a novel attempt to extract the seven ungiven parameters of PEMFCs stack. The sum of squared deviations (SSDs) among the measured and the relevant model-based calculated datasets is adopted to define the cost function. A Kepler Optimization Algorithm (KOA) is employed to decide the best values of these parameters within viable ranges. Initially, the KOA-based methodology is applied to assess the steady-state performance for four practical study cases under several operating conditions. The results of the KOA are appraised against four newly challenging algorithms and the other recently reported optimizers in the literature under fair comparisons, to prove its superiority. Particularly, the minimum values of the SSDs for Ballard Mark, BCS 0.5 kW, NedStack PS6, and Temasek 1 kW PEMFCs stacks are 0.810578 V2, 0.0116952 V2, 2.10847 V2, and 0.590467 V2, respectively. Furthermore, the performance measures are evaluated on various metrics. Lastly, a simplified trial to upgrade Amphlett’s model to include the PEMFCs’ electrical dynamic response is introduced. The KOA appears to be viable and may be extended in real-time conditions according to the presented scenarios (steady-state and transient conditions).


Mathematical formulation of PEMFCs' model
As formerly-stated, Amphlett's model is regarded as the most powerful method to simulate the dependency of the PEMFC's terminal voltage on the load current through various steady-state cases [22][23][24] .In order to examine the PEMFCs stack in operation, Mann and Amphlett's PEMFCs model makes a number of simplifying assumptions 18,40 .Here are a few of these presumptions: (i) The fuel and oxidant gas compositions are constant and precisely blended throughout the entire cell, (ii) The electrodes completely consume the fuel and oxidant, preventing any buildup inside the cell, (iii) There is no gas mixing between the anode and cathode compartments because the fuel and oxidant gases run parallel to the electrodes, and (iv) The porosity and thickness of the material do not change over time, (v) The fuel and oxidant gases' temperature, pressure, and humidity remain constant across the entire cell, and (vi) the electrolyte is a perfect ion conductor with no limits on electronic conductivity or mass transfer.
In spite of the aforementioned, the Mann's model can be extended to describe the transient behavior as well 24 .Amphlett supposes that the PEMFC's terminal voltage is subjected to three voltage decays which are the activation, the ohmic, and the concentration voltage drops, as described in Fig. 1 29 .The next few statements concisely illustrate the model, as thoroughly described in the state-of-the-art.The terminal voltage per single cell V cl ( V ), is given by (1) 25 .
where, E n is the Nernst open-circuit voltage per PEMFC in ( V ), V ac is the activation over-potential in ( V ) as the startup chemical reactions are relatively slow, V oh is the ohmic voltage loss due to the total resistance of the membrane and external leads in ( V ), and V cn is the concentration over-potential due to the high water content in the membrane at heavy loading in ( V).
As earlier-mentioned, to attain a high terminal voltage, N PEMFCs are connected in series creating a PEMFCs' stack whose terminal voltage is V st ( V ) and calculated by (2) 26 .
Generally, (2) supposed that all the PEMFCs are identical and behaves the same 6 .The E n can be calculated by (3) for T ≤ 373.15K 27 .
where, the partial pressures of hydrogen and oxygen are represented by P H and P O in ( atm ) and determined by ( 4) and ( 5), respectively 25 .T is the PEMFCs' stack operating temperature in ( K).
where, the load current and membrane useful area are denoted by I and A m in ( A ) and ( cm 2 ), respectively.
RH a and RH c express the relative humidities of the water steam at the anode and cathode regions, respectively.The anode and the cathode inlet pressures are indicated by P a and P c in ( atm ), respectively.The saturated vapor pressure is symbolized by P w in ( atm ), which can be given by ( 6) 29 .
Furthermore, V oh can be expressed by (9)   where, the resistances of the membrane and leads are denoted by R m and R l in ( �) , respectively.l points out the membrane thickness in ( cm ).The membrane specific resistivity is represented by ρ m in ( cm ) and determined by (10)  33 .
(1)  www.nature.com/scientificreports/It's worth affirming that refers to the water content in the membrane and its determination is a challenging issue because of its reliance on the cell-drawn current.However, a certain water content, at all possible operating conditions, is assumed unchangeable in this work.
where, δ is an empirical factor in ( V ) and J mx is the peak current density in ( Acm −2 ).Lastly, the stack output power (P st ) can be determined by (12) 36,41 .
At this moment, the reader can easily notice the seven unknown parameters are namely ( ) that are ungiven in the manufacturers' datasheets [31][32][33][34][35][36][37]39 . Thee parameters are optimally generated using the KOA-based methodology to achieve an accurate simulation of the PEMFCs under various operating scenarios.
where, n defines the number of the measured voltage-current dataset points.
Furthermore, the SSDs is susceptible to inequality bounds where each unknown parameter has its own lower and higher limits.It's noteworthy to indicate that the KOA preserves these limits while searching for optimal values.Utilizing the SSDs values, the main goal is to significantly fit the recorded terminal voltages to the relevant computed ones by the KOA-based methodology.

Kepler optimization algorithm
Basically, KOA is a physics-based optimizer that imitates the planets' motion according to Kepler's laws.Particularly, the sun and its planets (objects) moving around it in (fictitious) oval paths (orbits) are utilized to simulate the search space, which represents Kepler's first statement.Specifically, the planets in KOA (nominated solutions) exist at various positions from the sun (optimal solution) and at different times, hence, the exploration and exploitation concepts are effectively performed, as shown in Fig. 2.There are many factors affecting the planet's position from the optimal solution (the sun) such as the actual planet's position, the attraction force between the sun and the planet, and its revolving speed around the sun 39 .
Principally, KOA starts with stochastic initialization of the objects' numbers and positions, as represented by (15).where Y j i denotes the i th object (nominated solution) in the search area.The population of the solution nominees in the search area is represented by N p and the number of variables to be optimized is defined by d .lb j i and hb j i are the lower and higher boundaries of the j th decision parameter, respectively.
The initialization of the orbital deviation O d for each i th object is given by ( 16).where s is a stochastically normal generated number.The gravitational force G f i that attracts any plant Y i to the sun Y s is defined by the universal gravitational law which is function of the sun and the planet mass m s and m i , respectively, and length between the star and the planet R i , as given by (18).
where the global gravitational constant is symbolized by α and formulated by (19), and ǫ is a trivial number to avoid dividing by zero error.R i defines the normalized value of R i which is calculated by (20).The normalized values of m s and m i are represented by m s and m i and computed by ( 21)- (24), respectively.
where, t and t m are the actual iteration and the highest number of iterations.The initial value and a constant are denoted by α o and γ , respectively.
where s 4 is a stochastic number generated from 0 to 1. f k (t) refers to the fitness function of the kth object.
According to Kepler's law, the velocity of a planet is dependent on its position from the sun.In detail, if the planet is close to the sun, it will experience a strong gravitational force, and hence, it will accelerate its motion to avoid pulling into the sun and vice versa.So, the planet's velocity is described by (25).( 17) where a i (t) symbolizes a semi-main axis of the oval path of the i th planet and is given by (36).
Generally, most objects move around the sun in an anti-clockwise direction, nevertheless, some of them may revolve in a clockwise motion.KOA emulates this phenomenon by employing a flag F to avoid trapping in the local minima zones.More specifically, KOA use F to adjust the search flow so that the objects enhance their scanning ability in the search space.
The exploration phase is attained when the objects are away from the sun, referring that KOA efficiently explores the whole search space.Thence, the updated position of each planet away from the sun is described by (37).
On the other hand, if the objects are close to the star, KOA will concentrate on optimizing the exploitation phase.The switching between the two phases is done using an adaptive controlling factor h , which gradually alters as a function of the time, as revealed in (38).Thus, the new position of the planets based on this controlling strategy is formulated by (39).
where a 2 is a periodical regulating factor gradually decays from -1 to -2 for M cycles throughout the overall optimization task, as indicated in (40).
Finally, the optimum position of the objects and the sun is determined by (41).
It's clearly noticed that the KOA has only to parameters need to be manually set, which are N p and t m .As a result, lower computational burden and less independents runs are required to obtain the best performance of KOA.The overall steps of the KOA are presented, in detail, in Fig. 3 39 .

Numerical simulations and algorithm verification
To appraise the robustness and the efficacy of the proposed KOA, the performance of four commercial stacks, commonly studied in the state-of-art, are assessed under different steady-state circumstances.In addition, the unknown parameters' limits have been extracted from the literature and maintained unchanged over all the test cases, to guarantee fair judgment to all the competitive algorithms.Moreover, all of the numerical simulations are carried out on a computer with an Intel Core i7 CPU and 8 GB of RAM using the MATLAB platform, version R2022a.The operating system is Windows 10 Enterprise.
It's worth asserting that the values of the KOA controlling parameters N p and t m are 10 and 10,000, respectively.Furthermore, the best values of unspecified parameters are picked after executing KOA over 20 independent runs due to the high randomness of MOAs.At the end, the statistical performance of KOA is checked to prove its robustness.

Test cases' datasheets and the parameters' boundaries
The reader is invited to browse Table 1 for the technical specifications of the four earlier-announced test cases, which are extracted from 25,[27][28][29][30] .The relative humidity's of the vapor at anode and cathode RH a and RH c are maintained at 1.00 in all cases.Also, Table 1 (last three columns) reveals the minimum and maximum operating boundaries of the undefined parameters, as obtained from 34,36,37,42 . Initialize

KOA-based parameters' estimation outcomes
At this moment, KOA besides, two recent and two well-matured optimizers, are employed to efficiently determine the seven ungiven parameters of Amphlett's model.Actually, the driving training-based algorithm (DTBA) 43 and the election-based optimization algorithm (EBOA) 44 are the two new algorithms, while grey wolf algorithm (GWA) and particle swarm algorithm (PSA) represent the two benchmark algorithms.Practically, Tables 2, 3, 4 and 5 elucidate the KOA-based minimum SSD's values, for the four aforementioned test cases after 20 autonomous runs, compared to the four executed algorithms and the other recently-reported optimizers.Examples are improved artificial hummingbird algorithm (IAHA) 45 , honey badger algorithm (HBA) 46 , manta rays foraging algorithm (MRFA) 47 , pathfinder algorithm (PFA) 48 , neural network algorithm (NNA) 49 , and moth-flame algorithm (MFA) 50 .Besides, sparrow search algorithm (SSA) 51 , vortex search algorithm (VSA) 52 , modified monarch butterfly algorithm (MMBA) 53 , quasi oppositional bonobo algorithm (QOBA) 54 , modified farmland fertility algorithm (MFFA) 55 , marine predator algorithm (MPA) 56 , modified AEA (MAEA) 57 , satin bowerbird algorithm (SBA) 58 , and shark smell optimizer (SSO) 59 are also brought to comparison with the proposed KOA-based results.Over and above, the optimal values of the unknown parameters for each algorithm are also captured in the above-stated tables.Furthermore, the convergence trends of the applied optimizers for the four PEMFCs' stacks are revealed in Fig. 4a-d.The observation we may mention about KOA achieving excellent results for the studied applications but being hampered by a slight delay in convergence (required high iteration) is indeed an important limitation compared to those published in the literature.Convergence refers to the point at which an algorithm's performance stabilizes, indicating that it has learned the underlying patterns in the data.1).So, an unapplicable solution.

Algorithms
Parameters A delay in convergence can impact the efficiency and effectiveness of an algorithm, particularly in real-time or time-sensitive applications where prompt decision-making is crucial.However, it's worth noting that convergence speed can vary depending on the complexity of the problem, the size of the dataset, and the specific algorithm being used.
It may be affirmed from the convergence patterns that KOA outperforms the other candidates in terms of escaping from getting trapped in local minima while maintaining the fastest rate to reach the best SSD throughout 10,000 iterations.Accordingly, the polarization (V-I) curves of the actually recorded and the KOA, GWA, PSA, DTBA, and EBOA-based simulated datasets for the four commercial PEMFCs' stacks are captured in Fig. 5a-d.A closer look at Fig. 5a-d, it can be caught that the computed V-I curves, generated after injecting the optimal Table 3. KOA outcomes compared to other competitive optimizers for BCS 0.5 kW.Significant values are in bold.

Algorithms
Parameters www.nature.com/scientificreports/values of the parameters to the model, are consistent and well-fitted to the relevant experimental values.Furthermore, the minor SSD's values support the afore-said statements (see Tables 2, 3, 4, 5).Another form of KOA verification is the percentage terminal voltage error V %TE , which is determined to appraise how reliable and accurate the KOA is to fit the model-based calculated terminal voltages to the experimental ones, as described by ( 42) 36,46 .Figure 6a-d elucidates the changeability of V %TE along with the relevant Table 5. KOA outcomes compared to other competitive optimizers for Temasek 1 kW.Significant values are in bold.Now, after validating the proposed KOA-based methodology, it's time to utilize its outcomes to study the impact of the polarization losses on the PEMFCs' voltage profile.As previously discussed, the terminal voltage of the PEMFCs' stack is influenced by three voltage drops; activation, ohmic, and concentration losses 54,57 .In this regard, Fig. 7a-d reveal the alternation of every single polarization loss besides, the total losses V TL as functions of the stack drawn current, for all test cases.It can be concluded from Fig. 7a-d that when starting the PEMFCs at light load, the activation losses rapidly increase, then at intermediate loading values, it almost gets saturated, while the ohmic losses start a linear rise.At heavy drawn currents, the concentration losses considerably arise.

KOA-based outcomes under various operating conditions
In this subsection, after accrediting the KOA-based results, the effect of varying the PEMFCs' operating factors, P O /P H and T , on the polarization characteristics (V-I and P-I curves) can be adequately evaluated.Particularly, only two test cases are addressed for this purpose to avoid a lengthy article.More specifically, The P-I and V-I curves of BCS 0.5 kW and NedStack 6 kW under changing the temperature (40, 60, and 80 °C), while maintaining the other factors constant, are caught in Figs.8a,b and 9a,b, respectively.Moreover, Fig. 8c,d illustrates the impact of adjusting the suppliant partial pressures ( P O /P H = 0.2095/1, 1/1.5, 1.5/2.5 bar) on the aforementioned curves of BCS 0.5 kW, while unchanging the other datasheet' factors.The same is done for NedStack 6 kW, the V-I and P-I curves are generated while varying the input partial pressures ( P O /P H = 1/1, 1.5/2, 2/3 bar), as revealed in Fig. 9c,d 25,36 .
It goes without saying that according to Figs. 8 and 9, increasing the operating factors of the PEMFCs, whether P O /P H or T , positively affects the stack output voltage and power at the same drawn current.However, any increase over the nominal limits stated by the manufacturer will lead to severe degeneration of PEMFC's performance.for the activation losses.Thence, the mathematical representation of the dynamic activation voltage drop V ac (t) is given by (43) 23,57 .
where, the time constant of the activation over-potential reaction is symbolized by t d .It's worth stating that this simplified version of Amphlett's model doesn't consider the time constants of the balance-of-plant devices, like the hydrogen and oxygen compressors 24,[60][61][62] .Based on the above hypotheses, a simplified dynamic model is constructed via SIMULINK environment to simulate the PEMFC's electrical performance due to load variations, as shown in Fig. 10.
As a representative of the other test cases, NedStack 6 kW stack is employed to check the validity of the proposed dynamic model, where t d equals 1.2 s and the other datasheet's parameters are kept at their nominal values (see Table 1 (fourth column)).Specifically, a step variation of the stack drawn current, whose values are equal to 9 A, 171 A, 45 A at 0 s, 40 s, and 70 s, respectively, is applied to the proposed model, as cropped in Fig. 11a.The dynamic model simultaneously responds to this variation, as revealed in Fig. 11b,c.
More specifically, the stack terminal voltage values due to these changes are captured in Fig. 11b, which totally fit to the simulated steady-state values indicated in Fig. 5c.Furthermore, the behaviors of the internal polarization losses due to the afore variation are depicted in Fig. 11c.It's noteworthy to recall what is stated before, as the load current increases, the stack polarization losses also increase, and accordingly its terminal voltage reduces until (43)  V ac (t) = V ac (see( 7)) × 1

Conclusion and future prospective
A viable and effective methodology based on KOA has been introduced for optimally specifying the undefined parameters of the commonly investigated PEMFC's model, called Amphlett's model.Four practical study cases of well-known commercial PEMFCs' types have been comprehensively discussed through a set of simulated electrical characteristics such as the calculated V-I and P-I curves and the polarization losses alternation with the drawn currents.In addition to that it's worth mentioning that the percentage terminal voltage error between the KOA-based computed voltages and the experimentally recorded ones are 2.7%, 0.49%, − 1.33%, − 3.99% for Ballard Mark V, BCS 0.5 kW, NedStack PS6, and Temasek 1 kW PEMFCs, respectively.Furthermore, the effect of varying the input operating factors of the PEMFCs, like the temperature and the suppliant pressures, is deeply assessed for the whole test cases.Statistically, KOA robustness and preciseness have been validated through multiple indices like StD, min, max, mean, and computational CPU time, where it extremely outperforms the other challenging competitors.Moreover, the electrical dynamic performance of the PEMFCs is brought under study using an upgraded version of Amphlett's model.This work still needs to be expanded in order to accurately analyze the performance of PEMFCs stack due to actual disturbances and other operating parameters like the hydrogen and oxygen flow rates.Additionally, the electrical behavior of the PEMFCs stacks when operating in parallel is taken into account in the anticipated future work.Once again, the promising results of the KOA encourage the research community to extend this current to study the performance of such PEMFCs stacks with maximum power tracker in real conditions and to study their behaviors' when they are connected to microgrid.

Figure 4 .
Figure 4. SSD convergence patterns for the test cases.

Figure 5 .
Figure 5. V-I curves for the four test cases.

t d .e −t t d + 1 Figure 7 .
Figure 7. Plots of polarization losses for the test cases.

Figure 8 .
Figure 8. V-I and P-I curves of BCS 0.5 kW under various operating conditions.

Figure 9 .
Figure 9. V-I and P-I curves of NedStack 6 kW under various operating conditions. 32 Lastly, the orbital interval O t for each ith object is initialized by(17).

Table 1 .
Technical specs of the test cases and practical limits of the undefined parameters.

Table 2 .
KOA outcomes compared to other competitive optimizers for Ballard Mark V. Cpt computed.Significant values are in bold.*The published values violate the earlier-mentioned boundaries (see Table

Table 4 .
KOA outcomes compared to other competitive optimizers for NedStack PS6.Significant values are in bold.

Table 6 .
Statistical indices of KOA and others.Significant values are in bold.